#include <cmath>
using namespace std;
#include "model.h"
#include "../../ThirdPartyCode/Utilities.h"

//const double PI = 3.141592653589793;

// =================================================================
model::model(int im,int iSiteNum, double idt)
{
	int i;

	// Parameters
	m          = im;
	dt         = idt;
	SiteNum    = iSiteNum;

	K_alpha    = new double[m];
	K_beta     = new double[m];
	double dx_EW = PI/SiteNum;
	for(i=0;i<m;i++)
	{
		K_alpha[i] = 4/(2*SiteNum*dx_EW*dx_EW*dx_EW)*sin(i*PI/(2*SiteNum))*sin(i*PI/(2*SiteNum));
		K_beta[i] = K_alpha[i];
	}
	
	// Initialize states
	meanAlpha2 = new double[m];
	meanBeta2  = new double[m];
	varAlpha2  = new double[m];
	varBeta2   = new double[m];
	for(i=0;i<m;i++)
	{
		meanAlpha2[i] = 0.0;
		meanBeta2[i]  = 0.0;
		varAlpha2[i]  = 0.0;
		varBeta2[i]   = 0.0;
	}
	h          = 0.0;
	rho        = 0.0;
	time       = dt;
}

model::~model()
{
	delete [] meanAlpha2;
	delete [] meanBeta2;
	delete [] varAlpha2;
	delete [] varBeta2;
	delete [] K_alpha;
	delete [] K_beta;
}

model& model::operator=(const model& rhs)
{
	if (this == &rhs)
		return *this;

	int i;

	m          = rhs.m;
	dt         = rhs.dt;
	meanAlpha2 = new double[m];
	meanAlpha2 = new double[m];
	meanBeta2  = new double[m];
	varAlpha2  = new double[m];
	varBeta2   = new double[m];
	K_alpha    = new double[m];
	K_beta     = new double[m];
	
	for(i=0;i<m;i++)
	{
		meanAlpha2[i] = rhs.meanAlpha2[i];
		meanBeta2[i]  = rhs.meanBeta2[i];
		varAlpha2[i]  = rhs.varAlpha2[i];
		varBeta2[i]   = rhs.varBeta2[i];
		K_alpha[i]    = rhs.K_alpha[i];
		K_beta[i]     = rhs.K_beta[i];
	}
	return *this;
}

model::model(const model& rhs)
{
	int i;
	
	m          = rhs.m;
	dt         = rhs.dt;
	meanAlpha2 = new double[m];
	meanAlpha2 = new double[m];
	meanBeta2  = new double[m];
	varAlpha2  = new double[m];
	varBeta2   = new double[m];
	K_alpha    = new double[m];
	K_beta     = new double[m];
	
	for(i=0;i<m;i++)
	{
		meanAlpha2[i] = rhs.meanAlpha2[i];
		meanBeta2[i]  = rhs.meanBeta2[i];
		varAlpha2[i]  = rhs.varAlpha2[i];
		varBeta2[i]   = rhs.varBeta2[i];
		K_alpha[i]    = rhs.K_alpha[i];
		K_beta[i]     = rhs.K_beta[i];
	}
}
// =================================================================
void model::setState(struct lattice_state iState)
{
	int i;
	h = iState.H;
	rho = iState.SOR;
	for(i=0;i<m;i++)
	{
		meanAlpha2[i] = iState.alphas[i]*iState.alphas[i];
		meanBeta2[i]  = iState.betas[i]*iState.betas[i];
		varAlpha2[i]  = 0;
		varBeta2[i]   = 0;
	}
}
// =================================================================
void model::reset()
{
	int i;
	time = dt;
	for(i=0;i<m;++i)
	{
		meanAlpha2[i] = 0.0;
		meanBeta2[i]  = 0.0;
		varAlpha2[i]  = 0.0;
		varBeta2[i]   = 0.0;
	}
	h = 0.0;
	rho = 0.0;
}
// Update state variables ==========================================
void model::update(double sysTime,double w,double T)
{
	if(sysTime<time)
		return;
	else{
		int    j,n;
		double K;
		double tau;
		double nu;
		double sigma2;
		double rh;
		double temp1,temp2,temp3;
		double h_old,rho_old;
		// Save old state -----------------------------------
		h_old   = h;
		rho_old = rho;
		// k,tau,rh,nu,sigma2 -------------------------------
		K      = interp1(TAB_T,TAB_K,TAB_LENG,T);
		tau    = interp1(TAB_T,TAB_TAU,TAB_LENG,T);
		nu     = interp1(TAB_T,TAB_NU,TAB_LENG,T);
		sigma2 = interp1(TAB_T,TAB_SIGMA2,TAB_LENG,T);
		rh     = interp1(TAB_T,TAB_RH,TAB_LENG,T);
		// H ------------------------------------------------
		h = h_old+dt*rh;
		
		// rho ----------------------------------------------
		rho=(rho_old*h_old+rh*(K*dt+(K-rho_old)*tau*(exp(-dt/tau)-1)))/(h_old+dt*rh);
		
		// meanAlpha2, meanBeta2 ----------------------------
		for(j=0;j<m;j++)
		{
			n = j+1;
			temp1 = sigma2/(2*nu*n*n*dt);
			temp2 = exp(-2*nu*n*n*dt);
			meanAlpha2[j] = temp1+(meanAlpha2[j]-temp1)*temp2;
			meanBeta2[j]  = temp1+(meanBeta2[j]-temp1)*temp2;
		}
		
		// varAlpha2, varBeta2 ------------------------------
		for(j=0;j<m;j++)
		{
			n = j+1;
			temp1 = exp(-2*nu*n*n*dt);
			temp2 = exp(-4*nu*n*n*dt);
			temp3 = sigma2*(temp1-1)/(-2*nu*n*n);
			varAlpha2[j] = temp2*varAlpha2[j]+4*temp1*temp3*meanAlpha2[j]+2*temp3*temp3;
			varBeta2[j]  = temp2*varBeta2[j]+4*temp1*temp3*meanBeta2[j]+2*temp3*temp3;
		}
		time = time+dt;
	}
}
void model::update(double W)
{
	;
}
// Output rountines ============================================
double model::getR2()
{
	int j;
	double R2 = 0;
	for(j=0;j<m;j++)
	{
		R2 += (meanAlpha2[j]+meanBeta2[j])/(2*PI);
	}
	return R2;
}
// -------------------------------------------------------------
double model::getM2()
{
	int j;
	double M2 = 0;
	for(j=0;j<m;j++)
	{
		M2 += (meanAlpha2[j]*K_alpha[j]+meanBeta2[j]*K_beta[j])/(2*SiteNum);
	}
	return M2;
}
// -------------------------------------------------------------
double model::getVarR2()
{
	int j;
	double varR2 = 0;
	for(j=0;j<m;j++)
	{
		varR2 += (varAlpha2[j]+varBeta2[j])/(4*PI*PI);
	}
	return varR2;
}
// -------------------------------------------------------------
double model::getH()
{
	return h;
}
// -------------------------------------------------------------
double model::getRho()
{
	return rho;
}
// =============================================================
